Evaluate marine heatwaves(MHW)
cywhale
2022-06-16
R code modified from https://github.com/mjacox/Thermal_Displacement/blob/master/thermal_displacement.m
Note: “00_OISST_include_plotx.R” contains included library and gplotx function in previous Rmd files.
if (!require("sp")) install.packages("sp"); library(sp)
if (!require("sf")) install.packages("sf"); library(sf)
hdt <- read_stars("../data_src/oisst/monthly_heatwave_detrend/202002_heatwave.nc") %>% as.data.table() %>%
.[,`:=`(lng=fifelse(x>180, x-360, x), lat=y)] %>% .[,.(lng, lat)]
latx0 <- sort(unique(hdt$lat))
lngx0 <- sort(unique(hdt$lng))
Detrend <- TRUE ## selece result in monthly_heatwave_detrend
RE <- 6371 #km of Earth radius
res <- lngx0[2]-lngx0[1]
bbox= c(115.0, 20.0, 135.0, 35.0) #c(-180.0, -90.0, 180.0, 90.0) #first try a smaller bounding box
grd_x <- res
grd_y <- res
lngx1 <- seq(as.integer(-360*1000+grd_x*1000),
as.integer(360*1000-grd_x*1000), by=as.integer(grd_x*1000))/1000.0
latx1 <- rev(seq(as.integer(-90*1000+(grd_x/2)*1000),
as.integer(90*1000), by=as.integer(grd_x*1000))/1000.0)
res
lngxt <- sort(fifelse(lngx0<0, lngx0+360, lngx0))
latxt <- rev(latx0) #in NetCDF, lat is from 89.875 to -89.875
#lngxt
#latxt
dt <- fread("unzip -p ../data_src/mapfiles/sea_icemask_025d.csv.zip")
#dt
setorder(dt, -y, x)
dt[,ry:=rowid(x)] ## NOTE that in OISST.nc file, y is from Northest (89.875) to Southest (-89.875)
#################### Actually, x is 0-360 in nc so is 0-180 -180-0 arranged
mask <- dcast(dt, x ~ ry, value.var = "seamask")[,-1] %>% as.matrix()
mlon <- dcast(dt, x ~ ry, value.var = "lon")[,-1] %>% as.matrix()
mlat <- dcast(dt, x ~ ry, value.var = "y")[,-1] %>% as.matrix()
if (Detrend) {
indir <- "../data_src/oisst/monthly_heatwave_detrend/"
andir <- "../data_src/oisst/monthly_anom_detrend/"
} else {
indir <- "../data_src/oisst/monthly_heatwave/"
andir <- "../data_src/oisst/monthly_anom/"
}
td_max = 3000 # km
latn <- length(latx1) #dim(dm)[2] #720
dlon <- length(lngx1) #dim(dm)[3] #2879 (difference of longitude)
lonn <- dim(mask)[1] #1440
yrng <- seq(1982,2022)
endt <- "2022-05-22"
trackdate <- seq.Date(as.IDate(endt)-6, as.IDate(endt), by="day")
curryr <- year(as.IDate(endt))
currmo <- month(as.IDate(endt))
setkey(dt, x, y) #Note it would reorder dt to increasing y, and note that in dt, longitude is "lon" (old version is "lng")
rowGrp <- 3L
sf::sf_use_s2(FALSE) #Error in s2_geography_from_wkb(x).... in ppp_estimator(), for newer sf version
plan(multisession)
options(future.globals.maxSize= 1048576000*2)if (!require("raster")) install.packages("raster"); library(raster)
Rerasterize <- FALSE
if (Rerasterize) {
sp_coast <- ne_coastline() ## for ppp analysis, to check result
# for 0-360 degree, not break whole contour, so require a continuous world...
# https://gis.stackexchange.com/questions/266535/change-a-raster-from-longitude-display-180-180-to-0-360
sx <- raster::rasterize(sp_coast, raster::raster())
x1 <- crop(sx, extent(-180, 0, -90, 90))
x2 <- crop(sx, extent(0, 180, -90, 90))
x3 <- copy(x1)
extent(x1) <- c(180, 360, -90, 90)
mcoa <- merge(x3, x1, x2)
save(mcoa, file="../data_src/mapfiles/coast_extend.RData")
} else {
load(file="../data_src/mapfiles/coast_extend.RData")
}
plot(mcoa) #if want to plot in ppp_estimator, argument 'coast'=mcoa
if (!require("spatstat")) install.packages("spatstat"); library(spatstat)
if (!require("maptools")) install.packages("maptools"); library(maptools)
if (!require("rgeos")) install.packages("rgeos"); library(rgeos)
if (!require("ggrepel")) install.packages("ggrepel"); library(ggrepel)
if (!require("geojsonsf")) install.packages("geojsonsf"); library(geojsonsf)
## Null parameter in vector check
has_val <- function(x) {
if(is.null(x)) return(FALSE)
if(!length(x)) return(FALSE)
if(is(x, "data.frame")) return(nrow(x)>0) #"data.frame" %chin% class(x)
if(is.list(x)) x<-unlist(x, use.names = FALSE)
return(any(na.omit(x)!="") & any(!is.na(x)))
}
is_miss <- function(x) {!has_val(x)}
# code from odbapi::intersect_polyx/cywhale
intersect_poly <- function(site, poly, xy=c("longitude", "latitude"), polyid=NA, onlyInPoly=TRUE, crs=4326) {
sitex <- st_as_sf(site, coords = xy, crs = crs, remove = FALSE)
sjt <- st_intersects(sitex[,c(xy, "geometry")], poly)
setDT(sitex)[,overlap:=sapply(sjt, function(x) if(length(x)>1) {length(x)} else {1}, simplify = TRUE)][
,hx:=sapply(sjt, function(x) if(length(x)) {x} else {NA}, simplify = TRUE)
][,idtstx:=.I]
dt <- sitex[rep(1:.N, overlap)][, gid:=NA_integer_][ #### Handle overlapping polygons
!is.na(hx), gid:=as.integer(hx[[rleid(idtstx)]]), by=.(idtstx)
][,`:=`(hx=NULL, geometry=NULL, idtstx=NULL)]
if (onlyInPoly) dt <- dt[!is.na(gid),]
if (is_miss(polyid)) return(dt) #### no specify which poly id in poly, return index of rows in poly
xpoly <- polyid[1]
if (xpoly %chin% colnames(poly)) { #### if polyid is colnames in poly, return actural polyid by polyid[gid]
dt[!is.na(gid), (xpoly):=get(xpoly, poly)[gid]][,gid:=NULL]
return(dt)
}
if (any(grep("id", colnames(poly)))) { #### if polyid only a desired colnames, not in poly, try to return poly$id[gid]
idx <- colnames(poly)[grep("id", colnames(poly))[1]]
dt[!is.na(gid), (xpoly):=get(idx, poly)[gid]][,gid:=NULL]
return(dt)
}
setnames(dt, chmatch("gid", colnames(dt)), xpoly)
return(dt) #### othewise, return index of rows in poly, change the name "gid" to desired polyid
}
jet.colors <-
colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan",
"#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))
ppp_estimator <- function(yr, mon, tdf=NULL, heatwave=NULL, coast=NULL, xPlot=FALSE, xGjson=TRUE) {
#yr=2004; mon=4; tdf=NULL; heatwave=NULL; coast=mcoa; xPlot=TRUE; xGjson=TRUE; #for debugging
jt <- as.integer(mon)
mon <- fifelse(jt<10, paste0("0",jt), paste0(jt))
if (is.null(tdf)) {
tdf <- paste0("../data_src/oisst/monthly_td/", yr, mon, "_td.csv")
}
if (is.character(tdf)) {
if (!file.exists(tdf)) {
print(paste0("Fail to output Year-month: ", yr," - ", mon, " because input file: ", tdf, " not exist!!"))
return(NULL)
}
zt <- fread(tdf) %>% .[td_flag==TRUE, .(lng, lat, sst, anom, xd, yd, sst_dis, ddis)]
} else if (is.data.frame(tdf)) {
zt <- setDT(tdf)[td_flag==TRUE, .(lng, lat, sst, anom, xd, yd, sst_dis, ddis)]
} else {
print(paste0("Fail to output Year-month: ", yr," - ", mon, " because input file or data not exist!!"))
return(NULL)
}
## Heatwaves occurs at zc with its sst, but copy a lng<0 to 180-360, make it extend from -180 - 0 - 180 - 360
zt[,`:=`(rowid=.I)]
zcm <- zt[lng<0, .(lng, lat, sst, rowid)]
zc1 <- zt[,.(lng, lat, sst, anom, rowid)]
zc <- st_as_sf(rbindlist(list(zc1[,.(lng,lat,sst,rowid)] %>% .[,lng:=fifelse(lng<0, lng+360, lng)], zcm), use.names = T), coords=c("lng", "lat"))
#st_crs(zc) <- 4326
ptc<- as.ppp(zc)
Window(ptc) <- owin(xrange=c(-180,360), yrange=c(-90,90))
kzc <- density(ptc, sigma=15, kernel = "quartic", adjust = 0.2, cut=6) #, ext=extent(x))
kzct <- copy(kzc)
kzct[kzct<5] <- NA
## Thermal displacement candidates occurs at zh with its sst (sst_dis)
zhm <- zt[xd<0, .(xd, yd, sst_dis, rowid)] %>% setnames(1:3, c("lng","lat","sst"))
zh1 <- zt[,.(xd, yd, sst_dis, ddis,rowid)] %>% setnames(1:3, c("lng","lat","sst"))
zh <- st_as_sf(rbindlist(list(zh1[,.(lng,lat,sst,rowid)] %>% .[,lng:=fifelse(lng<0, lng+360, lng)], zhm), use.names = T), coords=c("lng", "lat"))
pth<- as.ppp(zh)
Window(pth) <- owin(xrange=c(-180,360), yrange=c(-90,90))
kzh <- density(pth, sigma=15, kernel = "quartic", adjust = 0.2, cut=6)
kzht <- copy(kzh)
kzht[kzht<5] <- NA
if (xPlot) {
plot(kzc, main=NULL, las=1, box=FALSE,axes=FALSE, ribbon =FALSE, useRaster=FALSE, xlim=c(-180,360), ylim=c(-90,90))
contour(kzct, nlevels = 6, col="red", add=T)
contour(kzht, nlevels = 6, col="green",add=T)
if (!is.null(coast)) plot(coast, col="grey", add=T, legend=F, axes=F, box=F)
abline(h=0, col="lightgrey", lty=2, lwd=1)
}
#--------------------------------------------------------------------
## thermal displacement, need extend -180 - 180 - 360 to do raster analysis
gzc <- as(kzct, "SpatialGridDataFrame") # convert to spatial grid class
igzc<- as.image.SpatialGridDataFrame(gzc) # convert again to an image
cgzc <- contourLines(igzc, nlevels = 6) # create contour object - change 8 for more/fewer levels
sldc <- ContourLines2SLDF(cgzc, #proj4string= CRS("+init=epsg:4326")) #got the same as CRS(+proj=...)
CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")) # convert to SpatialLinesDataFrame
gzh <- as(kzht, "SpatialGridDataFrame") # convert to spatial grid class
igzh<- as.image.SpatialGridDataFrame(gzh) # convert again to an image
cgzh <- contourLines(igzh, nlevels = 6) # create contour object - change 8 for more/fewer levels
sldh <- ContourLines2SLDF(cgzh, CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")) # convert to SpatialLinesDataFrame
Polyc <- rgeos::gPolygonize(sldc[-1, ])
gac <- gArea(Polyc, byid = T)/10000
Polyc <- SpatialPolygonsDataFrame(Polyc, data = data.frame(gac), match.ID = F)
sf1 <- st_as_sf(Polyc) %>% st_union()
Polyc2 <- gPolygonize(sldh)
gah <- gArea(Polyc2, byid = T)/10000
Polyc2 <- SpatialPolygonsDataFrame(Polyc2, data = data.frame(gah), match.ID = F)
sf2 <- st_as_sf(Polyc2) %>% st_union()
if (xPlot) {
plot(sldc, col = heat.colors(8), main=NULL, las=1, axes=FALSE, xlim=c(-180,360), ylim=c(-90,90)) #add=T
plot(sldh, col = terrain.colors(8), add=T)
plot(sf1, col="#F8766D80", main=NULL, las=1, axes=FALSE, xlim=c(-180,360), ylim=c(-90,90), add=T)
plot(sf2, col="#00BFC480", add=T) #Polyc2
#plot(cen1, col="red", pch=10, add=T)
#plot(cen2, col="green", pch=7, add=T)
#abline(h=0, col="lightgrey", lty=2, lwd=1)
}
sf1x <- st_cast(sf1, "POLYGON") %>% st_sf(geom = .)
sf2x <- st_cast(sf2, "POLYGON") %>% st_sf(geom = .)
## Use sf, geojson output
xrng <- c(-90, 270) #c(0.360)
sflct <- st_as_sf(sldc, dim = "XY", crs=4326) %>% st_crop(xmin=xrng[1], ymin=-90, xmax=xrng[2], ymax=90) ##%>% ## clip [0, 360] range
# sflc$geometry <- (st_geometry(sflc) + c(360,90)) %% c(360) - c(0,90)
#sf::st_wrap_dateline(options = c("WRAPDATELINE=YES", "DATELINEOFFSET=180"))
sflc <- lapply(1:nrow(sflct), function(x) {
sf::st_wrap_dateline(sflct[x,],options = c("WRAPDATELINE=YES", "DATELINEOFFSET=180"))
}) %>% do.call(rbind, .)
st_crs(sflc) <- 4326
sflht <- st_as_sf(sldh, crs=4326) %>% st_crop(xmin=xrng[1], ymin=-90, xmax=xrng[2], ymax=90) #%>%
#sf::st_wrap_dateline(options = c("WRAPDATELINE=YES", "DATELINEOFFSET=180"))
sflh <- lapply(1:nrow(sflht), function(x) {
sf::st_wrap_dateline(sflht[x,],options = c("WRAPDATELINE=YES", "DATELINEOFFSET=180"))
}) %>% do.call(rbind, .)
st_crs(sflh) <- 4326
sf1y <- sf1x %>% #st_crop(xmin=xrng[1], ymin=-90, xmax=xrng[2], ymax=90) %>%
sf::st_wrap_dateline(options = c("WRAPDATELINE=YES", "DATELINEOFFSET=180"))
sf2y <- sf2x %>% #st_crop(xmin=xrng[1], ymin=-90, xmax=xrng[2], ymax=90) %>%
sf::st_wrap_dateline(options = c("WRAPDATELINE=YES", "DATELINEOFFSET=180"))
cen1 <- st_centroid(sf1x) #dim(sf1x) is n x 1 #Note that use continous sf1x, not sf1y(who take cake dateline)
cen2 <- st_centroid(sf2x) #dim(sf2x) is m x 1
if (FALSE) { #only Debug needed
xrng = range(c(st_coordinates(cen1)[,1], st_coordinates(cen2)[,1]))
yrng = range(c(st_coordinates(cen1)[,2], st_coordinates(cen2)[,2]))
plot(coast, col="grey", legend=F, axes=F, box=F, xlim=c(floor(xrng[1]), floor(xrng[2]+1)), ylim=c(floor(yrng[1])-1, floor(yrng[2]+1)+1)) #add=T,
plot(sf1x, border="#F8766D80", xlim=c(floor(xrng[1]), floor(xrng[2]+1)), ylim=c(floor(yrng[1]), floor(yrng[2]+1)), add=T)
text(st_coordinates(cen1), labels=seq_along(1:nrow(cen1)), col="red", cex=0.75)
plot(sf2x, border="#00BFC480", add=T)
text(st_coordinates(cen2), labels=seq_along(1:nrow(cen2)), col="blue", cex=0.75, adj = c(0, 1))
}
# st_nearest_feature auto conver long>180 to neg, and can be diff with nrx-> nmin result but still resonable
# nrx <- st_nearest_feature(cen1, cen2) ## cannot solve bug GEOS has differ version in Linking (3.7.1) & runtime (3.4.2)
nrx <- st_distance(cen1, cen2) #dim(nrx) is nxm (distance matrix)
nmin<- apply(nrx, 1, which.min) #BY ROW, length(nmin) = n; means which cen2 of sf2x (sst_displacement, DIS) has min distance to each cen1 (heatwaves, HW) of sf1x
mmin0 <- apply(nrx, 2, which.min) #BY Column, means which cen1 (heatwaves, HW) has min distance to each cen2 (displacement, DIS)
xmin<- mmin0
xmin[nmin]<- seq_along(nmin) # if cen1(n) -> cen2(m) can multiple-to-1 but should not 1-to-multiple
# Basically, we want to find which DIS has min dist to HW, multple HW can assign to the same DIS.
############################################################################################################
#cen1= #nmin= | #cen2 #mmin0 #xmin[nmin] xmin is a new estimator for cen1,
#seq_along( #which.min | #1 #25 #xmin[2] <- 1
#nmin) # of m cell | #2 #1 #xmin[22] <- 2
#1 # 2 | #.....................
#2 # 22 | #............#xmin[47] <- n
#.......................| #m #23
#n # 47 |
############# e.g min distance of (HW, DIS) is (1, 2) then we replace xmin[2] <- 1, so (DIS, HW) is (2, 1)
############# e.g min distance of (HW, DIS) is (2, 22) then we replace xmin[22]<- 2, so (DIS, HW) is (22,2)
# other not replaced in xmin is originally filled by mmin0, so we got some (DIS, with some duplicated HW)
# so unique(xmin[which(duplicated(xmin))] get those duplicated HW by unique <- smaller duplicated index, that is filled by seq_alon(nmin)
# then fetch value by nmin[some of seq_alon(nmin) caused from duplicated HW] === which.min of DIS cells to HW, as DUP2
# then unqiue(sort(NOT_DUP, DUP2))
if (any(duplicated(xmin))) {
munq <- unique(sort(c(which(!duplicated(xmin)), nmin[unique(xmin[which(duplicated(xmin))])]))) # select unique xmin, Note: xmin content is n, but munq is index of xmin
} else {
munq <- seq_along(xmin)
} # still will have duplicated(xmin[munq])
#zht <- copy(zh1)[,lng:=fifelse(lng<0, lng+360, lng)] #Nope, sf2x extend from -180 to 360 not only (0,360)
mply <- intersect_poly(zh1, sf2y[munq,], xy=c("lng", "lat"), polyid=NA, onlyInPoly=TRUE, crs=4326)
mply[,gid:=munq[gid]][,overlap:=NULL]
mmed <- mply[!is.na(gid), .(dmed=median(ddis)), by=.(gid)][,`:=`(nid=xmin[gid])]
jy <- xmin[munq]
if (any(duplicated(jy))) {
jm <- c(which(duplicated(jy)), which(duplicated(jy, fromLast = TRUE))) #index to xmin[munq] which is duplicated in its content (n)
selm <- mmed[nid %in% jy[jm], .(selm=gid[which.min(dmed)]), by=.(nid)]$selm
mmed <- mmed[gid %in% sort(c(munq[seq_along(jy)[-jm]], selm)),]
}
nunq <- sort(unique(mmed$nid))
nply <- #odbapi::intersect_polyx
intersect_poly(zc1, sf1y[nunq,], xy=c("lng", "lat"), polyid=NA, onlyInPoly=TRUE, crs=4326)
nply[,`:=`(nid=nunq[gid])][,`:=`(overlap=NULL, gid=NULL)]
setnames(mply, 1:3, c("xd","yd","sst_dis"))
mnx <- merge(mmed, merge(mply, nply, by=c("rowid"), all.x=T) %>% .[!is.na(nid),], # & sst>sst_dis,] <- MUST BE
by=c("gid","nid"), all=T) %>% .[!is.na(dmed),] %>% .[ ## just find an example in each (gid, nid)
,.(lng=lng[which.max(ddis)], lat=lat[which.max(ddis)],
sst=sst[which.max(ddis)], anom=anom[which.max(ddis)],
xd=xd[which.max(ddis)], yd=yd[which.max(ddis)],
sst_dis=sst_dis[which.max(ddis)], ddis=ddis[which.max(ddis)]), by=.(gid, nid, dmed)]
mnx[,`:=`(xm1=fifelse(sign(lng)*sign(xd) == -1 & (abs(lng)>150 | abs(xd)>150), fifelse(sign(lng) == -1, -180.0, 180.0), NA_real_), ym1=0.5*(lat+yd))][
,`:=`(xm2=fifelse(sign(lng)*sign(xd) == -1 & (abs(lng)>150 | abs(xd)>150), fifelse(sign(lng) == -1, 180.0, -180.0), NA_real_), ym2=ym1)]
sfl1 <- copy(sflc)
#sfl$level <- rev(as.integer(as.character(sflc$level))) - min(as.integer(as.character(sflh$level)))
sfl1$level <- c(1:nrow(sflc))
sfl1 <- st_zm(sfl1)
sfl2 <- copy(sflh)
#sfl2$level <- as.integer(as.character(sfl2$level))
sfl2$level <- c(1:nrow(sflh)) + 16 - nrow(sflh) #16 or 11 is the number or color palette allowed
sfl2 <- st_zm(sfl2)
sfl <- rbind(sfl1, sfl2)
if (xPlot) {
if (is.null(heatwave)) {
heatwave <- read_stars(paste0("../data_src/oisst/monthly_heatwave_detrend/", yr, mon, "_heatwave.nc"))
names(heatwave)[1] <- "heatwave"
}
ht <- as.data.table(heatwave) %>%
.[,`:=`(longitude=fifelse(x>180, x-360, x), latitude=y)] %>% .[,.(longitude, latitude, heatwave)]
g1 <- ggplot() +
geom_point(data=ht[heatwave==1,], aes_string(x="longitude", y="latitude"), alpha=0.65, col="gold", size = 0.01, stroke = 0, shape = 16) +
geom_point(data=zc1, aes(x=lng, y=lat), alpha=0.75, col="orange", size = 0.01, stroke = 0, shape = 16) +
geom_point(data=zh1, aes(x=lng, y=lat), alpha=0.75, col="purple", size = 0.01, stroke = 0, shape = 16) +
geom_sf(data = ne_coastline(scale = "large", returnclass = "sf"), color = 'darkgray', size = .3) +
geom_sf(data = sfl, aes(color=level)) + #factor(level)
#scale_colour_brewer(palette = "RdYlBu") +
#scale_colour_brewer(palette = "RdYlBu", direction = -1) +
scale_colour_gradientn(colours=rev(jet.colors(16))) +
geom_sf(data = sf1y, fill="#F8766D80", color="#F8766D80", alpha=0.5) +
geom_sf(data = sf2y, fill="#00BFC480", color="#00BFC480", alpha=0.5) +
guides(fill="none", color="none") + coord_sf() + #xlim(c(-180, 180)) + ylim(c(-90, 90)) +
scale_x_continuous(limits = c(-180, 180), expand = c(0, 0)) +
scale_y_continuous(limits = c(-90, 90), expand = c(0, 0)) +
labs(x="Longitude",y="Latitude") +
theme(
panel.background = element_rect(fill="white"),
panel.border = element_rect(colour = "black", fill=NA, size=0.75),
panel.grid.major = element_line(colour = "lightgrey"),
panel.grid.minor = element_line(colour = "lightgrey"),
strip.background = element_blank(),
strip.text.y = element_text(angle = 0),#,face = "italic"),
legend.background = element_rect(fill = "transparent", colour = "transparent"),#, #"white"),
legend.position = NULL
)
if (any(!is.na(mnx$xm1))) {
g1 <- g1 +
geom_segment(data=mnx[!is.na(xm1),], aes(x=lng, xend=xm1, y=lat, yend=ym1), size = 0.75, color="purple") + #, color=factor(seamask)
geom_segment(data=mnx[!is.na(xm1),], aes(x=xm2, xend=xd, y=ym2, yend=yd), size = 0.75, arrow = arrow(length = unit(0.08, "cm")), color="purple") +
geom_segment(data=mnx[is.na(xm1),], aes(x=lng, xend=xd, y=lat, yend=yd), size = 0.75, arrow = arrow(length = unit(0.08, "cm")), color="purple") +
geom_text_repel(data=mnx[!is.na(xm1),], aes(x=xd, y=yd, label=paste0(round(ddis,0), " km")), color = "black", size=3, hjust=1, vjust=1) + #geom_text
geom_text_repel(data=mnx[is.na(xm1),], aes(x=0.5*(lng+xd), y=0.5*(lat+yd), label=paste0(round(ddis,0), " km")), color = "black", size=3, hjust=1, vjust=1)
} else {
g1 <- g1 +
geom_segment(data=mnx, aes(x=lng, xend=xd, y=lat, yend=yd), size = 0.75, arrow = arrow(length = unit(0.08, "cm")), color="purple") +
geom_text_repel(data=mnx, aes(x=0.5*(lng+xd), y=0.5*(lat+yd), label=paste0(round(ddis,0), " km")), color = "black", size=3, hjust=1, vjust=1)
}
#dev.off()
#g1
}
# ==== Output as geojson
if (xGjson) {
outdir <- "../data_src/oisst/thermal_displace/"
sfc1 <- st_as_sf(zc1[,.(sst, anom, lng, lat)], coords=c("lng", "lat"), crs=4326) ## add sst anom too big for geojson 15Mb
jfc1 <- sf_geojson(sfc1)
geojson_write(jfc1, file = paste0(outdir, yr, mon, "_hwv_occ.geojson"))
sfc2 <- st_as_sf(zh1[,.(sst, ddis, lng, lat)], coords=c("lng", "lat"), crs=4326) ## add sst anom too big for geojson 15Mb
jfc2 <- sf_geojson(sfc2)
geojson_write(jfc2, file = paste0(outdir, yr, mon, "_hwv_dis.geojson"))
geojson_write(sf_geojson(sflc), file = paste0(outdir, yr, mon, "_hwv_occ_contour.geojson"))
geojson_write(sf_geojson(sflh), file = paste0(outdir, yr, mon, "_hwv_dis_contour.geojson"))
sf1d <- st_as_sf(data.frame(ID=seq_along(sf1y), geometry=st_geometry(sf1y)), crs=4326)
sf2d <- st_as_sf(data.frame(ID=seq_along(sf2y), geometry=st_geometry(sf2y)), crs=4326)
geojson_write(sf1d, file = paste0(outdir, yr, mon, "_hwv_occ_poly.geojson"))
geojson_write(sf2d, file = paste0(outdir, yr, mon, "_hwv_dis_poly.geojson"))
mnxf <- st_as_sf(mnx, coords=c("lng", "lat"), crs=4326)
geojson_write(mnxf, file = paste0(outdir, yr, mon, "_hwv_dis_correspond.geojson"))
}
print(paste0("Successfully output Year-month: ", yr," - ", mon, " geojson for thermal displacement"))
if (xPlot) {
return(g1)
}
}
apply_oisst_masks <- function(ii, jj, d_mask) { #, mask, mlat, mlon) {
# ====================================================
# Inputs:
# ii: longitude index on OISST grid
# jj: latitude index on OISST grid
# d: matrix of distances to all other OISST grid cells
# from point [ii, jj]. Dimensions are [lon lat]
# mask: mask defining regions, created by make_oisst_masks.m
# Dimensions are [lon lat]
# lat: Matrix of OISST latitudes
# lon: Matrix of OISST longitudes
# Output:
# d_mask: matrix of distances with unavailable locations masked out
# d_mask <- d
# Exclude ice-surrounded areas
d_mask[mask==3] <- NA_real_
if (!is.na(mask[ii, jj])) {
maskij <- letters[as.integer(mask[ii, jj])]
# Handle regional cases
switch(maskij,
d = { d_mask[mask != 4] <- NA_real_ }, # maskij == 4
e = { d_mask[mask != 5] <- NA_real_ }, # maskij == 5
f = {
d_mask[mask<=5 | mask==7 | mask==8 | mlat>48 | (mlat>43 & mlon>351)] <- NA_real_
if (mlon[ii,jj]>12 & mlon[ii,jj]<20 & mlat[ii,jj]>42.3 & mlat[ii,jj]<46) {
d_mask[mask!=6] <- NA_real_
}
}, # maskij == 6
g = { d_mask[!(mask==7 | mask==11)] <- NA_real_ }, # maskij == 7
h = { d_mask[!(mask==8 | mask==9)] <- NA_real_ }, # maskij == 8
i = { d_mask[!(mask==7 | mask==8 | mask==9 | mask==11)] <- NA_real_ }, # maskij == 9
j = { d_mask[!(mask==10 | mask==11)] <- NA_real_ }, # maskij == 10
k = { d_mask[(mask>=4 & mask<=6) | mask==12] <- NA_real_ }, # maskij == 11
l = { d_mask[mask>=9 & mask<=11] <- NA_real_ }, # maskij == 12
m = { d_mask[!(mask==13 | mask==14)] <- NA_real_ }, # maskij == 13
n = { d_mask[(mask>=15 & mask<=17) | mlat>283] <- NA_real_ }, # maskij == 14
o = { d_mask[mask==2 | mask==13 | mask==14 | mask==17 | mlon<260 | mlon>280] <- NA_real_ }, # maskij == 15
p = { d_mask[mask==13 | mask==14 | mlon<260] <- NA_real_ }, #maskij == 16
q = { d_mask[mask==15] <- NA_real_ } # maskij == 17
)
}
return(d_mask)
}Recalculate_displace <- FALSE
Rewrite_displace <- FALSE
Reestimate <- FALSE
if (Recalculate_displace) {
for (i in yrng) {
mmx <- fifelse(i==curryr, currmo-1L, 12L)
for (j in 1:mmx) {
monj <- fifelse(j<10, paste0("0",j), paste0(j))
print(paste0("Now in Year-month: ", i," - ", monj, " to calculate thermal displacement"))
filet <- paste0("../data_src/oisst/monthly_td/", i, monj, "_td.csv")
if (Rewrite_displace | (!file.exists(filet))) {
z <- read_stars(paste0(indir, i, monj, "_heatwave.nc"))
names(z)[1] <- "heatwave"
sst <- read_stars(paste0("../data_src/oisst/monthly_sst/", i, monj, "_sst.nc"))
names(sst)[1] <- "sst"
anom <- read_stars(paste0(andir, i, monj, "_anom.nc"))
names(anom)[1] <- "anom"
zd <- merge(as.data.table(z) %>% setkey(x,y), dt[,.(x,y,lon,seamask)], by=c("x","y"), all=T)[
heatwave==1 & (is.na(seamask) | seamask >3), .(lon, y, heatwave, seamask)][
,`:=`(xd=NA_real_, yd=NA_real_, ddis=NA_real_,
sst=NA_real_, anom=NA_real_, sst_dis=NA_real_, td_flag=NA,
rowgrp=cut(seq_len(.N), rowGrp, labels = FALSE))] #%>%
#.[sample.int(nrow(.), 3000),] ## Just test performance bottelneck
setnames(zd, 1:2, c("lng","lat")) ## Note here change "lon" -> "lng"
#print("Start get_tdx...")
get_tdx <- function(lng, lat) {
lngt <- fifelse(lng<0, lng+360, lng)
ii <- which(lngxt==lngt)
it <- which(lngx0==lng)
jt <- which(latx0==lat)
jj <- latn-jt+1 ## NetCDF lat in rev order
#if (is.na(mask[ii, jj]) | mask[ii, jj] > 3) { #z[[1]][ii, jj]==1 & () ## already check in zd[]
## (old version) Note that dm (aboving), diff of lng order is reversed to match NetCDF, not the same of ref matlab
## (old version use whole dm 11GB too large)
## d_lat <- apply_oisst_masks(ii, jj, t(dm[jj, ,seq(it, it+lonn-1)])) #, mask, latm, lngm) #t(dm[jj, ,seq(dlon-lonn-it+2, dlon-it+1)])
rlng <- seq(as.integer((0.125-lngt)*1000),
as.integer((359.875-lngt)*1000), by=as.integer(grd_x*1000))/1000.0
#rlng <- lngx1[seq(dlon-lonn-it+2, dlon-it+1)] #seq(it, it+lonn-1)
#dm <- array(rep(NA_real_, length(latx1)*length(rlng)), dim=c(length(latx1),length(rlng)))
dm <- t(Re(RE*acos(sin(latx1[jj]*pi/180)*sin(latx1*pi/180) +
cos(latx1[jj]*pi/180)*cos(latx1*pi/180) %*% t(cos(rlng*pi/180)))))
# print(object.size(dm), units="Mb") ## 7.9 Mb << 11Gb !!
# print(paste0("After matrix multiplication: ", rowid))
d_lat <- apply_oisst_masks(ii, jj, dm) #, mask, latm, lngm) #t(dm[jj, ,seq(dlon-lonn-it+2, dlon-it+1)])
sst_norm <- sst[[1]][ii,jj] - anom[[1]][ii,jj]; # i.e. climatology, because sst - climatology = anomaly
d_lat[is.na(sst[[1]]) | sst[[1]]>sst_norm | d_lat>td_max] <- NA_real_;
idx <- which(!is.na(d_lat))
if (any(idx)) {
minidx <- which.min(d_lat[idx])
xt <- idx[minidx] %% lonn
yt <- as.integer(idx[minidx]/lonn)
yt <- fifelse(xt==0, yt, yt+1)
xt <- fifelse(xt==0, lonn, xt)
xx <- lngxt[xt]
xx <- fifelse(xx>180, xx-360, xx)
yy <- latxt[yt]
return(list(xd=xx, yd=yy, ddis=d_lat[idx[minidx]],
sst=sst[[1]][ii,jj], anom=anom[[1]][ii,jj], sst_dis=sst[[1]][xt,yt], td_flag=TRUE)) #]
} else {
return(list(xd=NA_real_, yd=NA_real_, ddis=NA_real_,
sst=sst[[1]][ii,jj], anom=anom[[1]][ii,jj], sst_dis=NA_real_, td_flag=FALSE)) #]
}
#}
}
#for (ii in seq_along(lngx0)) {
# for (jj in seq_along(latx0)) {
zdx <- rbindlist(future_lapply(seq_len(rowGrp), function(grp) {
x <- zd[rowgrp==grp, ]
return(
x[,c("xd", "yd", "ddis", "sst", "anom", "sst_dis", "td_flag"):=get_tdx(lng, lat), by = 1:nrow(x)]
)
}), use.names = TRUE, fill = TRUE)
#} #use rowGrp =3 with future_lapply for nrow=3000 user system elapsed 0.742 0.033 70.476
#use rowGrp =4 for nrow=3000 user system elapsed 0.867 0.048 60.149
#For one month data 1400x720 rows need user system elapsed 21.244 1.086 3309.989
#zd[#heatwave==1 & (is.na(seamask) | seamask >3)
# , c("xd", "yd", "ddis", "sst", "anom", "sst_dis", "td_flag"):=get_tdx(lng, lat, rowid), by = 1:nrow(zd)] #by=.(lng, lat)]
setkey(zdx, lng, lat)
fwrite(zdx[,rowgrp:=NULL], file= filet)
#} #for-loop is very slow
}
if (Reestimate) {
ppp_estimator(yr=i, mon=j, xPlot=FALSE, xGjson=TRUE)
}
}
}
}g1 <- ppp_estimator(yr=2020, mon=2, xPlot=TRUE, xGjson=FALSE)
g1
yr=2018
mon="06"
outdir <- "../data_src/oisst/thermal_displace/"
sf1d <- st_read(paste0(outdir, yr, mon, "_hwv_occ_poly.geojson"))
sf2d <- st_read(paste0(outdir, yr, mon, "_hwv_dis_poly.geojson"))
sfl1 <- st_read(paste0(outdir, yr, mon, "_hwv_occ_contour.geojson"))
sfl2 <- st_read(paste0(outdir, yr, mon, "_hwv_dis_contour.geojson"))
sfc1 <- st_read(paste0(outdir, yr, mon, "_hwv_occ.geojson"))
sfc2 <- st_read(paste0(outdir, yr, mon, "_hwv_dis.geojson"))
mnxf <- st_read(paste0(outdir, yr, mon, "_hwv_dis_correspond.geojson"))
mnx <- cbind(as.data.table(sf::st_coordinates(mnxf)) %>% setnames(1:2, c("lng","lat")), as.data.table(mnxf))
mnx[xm1=="NA", xm1:=NA]
mnx[xm2=="NA", xm2:=NA]
sfl1$level <- c(1:nrow(sfl1))
sfl2$level <- c(1:nrow(sfl2)) + 16 - nrow(sfl2)
sflt <- rbind(sfl1, sfl2)heatwave <- read_stars(paste0("../data_src/oisst/monthly_heatwave_detrend/", yr, mon, "_heatwave.nc"))
names(heatwave)[1] <- "heatwave"
ht <- as.data.table(heatwave) %>%
.[,`:=`(longitude=fifelse(x>180, x-360, x), latitude=y)] %>% .[,.(longitude, latitude, heatwave)]
gf <- ggplot() +
geom_point(data=ht[heatwave==1,], aes_string(x="longitude", y="latitude"), alpha=0.65, col="gold", size = 0.01, stroke = 0, shape = 16) +
geom_sf(data=sfc1, alpha=0.75, col="orange", size = 0.01, stroke = 0, shape = 16) +
geom_sf(data=sfc2, alpha=0.75, col="purple", size = 0.01, stroke = 0, shape = 16) +
geom_sf(data = ne_coastline(scale = "large", returnclass = "sf"), color = 'darkgray', size = .3) +
geom_sf(data = sflt, aes(color=level)) +
scale_colour_gradientn(colours=rev(jet.colors(16))) +
geom_sf(data = sf1d, fill="#F8766D80", color="#F8766D80", alpha=0.5) +
geom_sf(data = sf2d, fill="#00BFC480", color="#00BFC480", alpha=0.5) +
guides(fill="none", color="none") + coord_sf() + #xlim(c(-180, 180)) + ylim(c(-90, 90)) +
scale_x_continuous(limits = c(-180, 180), expand = c(0, 0)) +
scale_y_continuous(limits = c(-90, 90), expand = c(0, 0)) +
labs(x="Longitude",y="Latitude") +
theme(
panel.background = element_rect(fill="white"),
panel.border = element_rect(colour = "black", fill=NA, size=0.75),
panel.grid.major = element_line(colour = "lightgrey"),
panel.grid.minor = element_line(colour = "lightgrey"),
strip.background = element_blank(),
strip.text.y = element_text(angle = 0),#,face = "italic"),
legend.background = element_rect(fill = "transparent", colour = "transparent"),#, #"white"),
legend.position = NULL
)
if (any(!is.na(mnx$xm1))) {
gf <- gf +
geom_segment(data=mnx[!is.na(xm1),], aes(x=lng, xend=xm1, y=lat, yend=ym1), size = 0.75, color="purple") + #, color=factor(seamask)
geom_segment(data=mnx[!is.na(xm1),], aes(x=xm2, xend=xd, y=ym2, yend=yd), size = 0.75, arrow = arrow(length = unit(0.08, "cm")), color="purple") +
geom_segment(data=mnx[is.na(xm1),], aes(x=lng, xend=xd, y=lat, yend=yd), size = 0.75, arrow = arrow(length = unit(0.08, "cm")), color="purple") +
geom_text_repel(data=mnx[!is.na(xm1),], aes(x=xd, y=yd, label=paste0(round(ddis,0), " km")), color = "black", size=3, hjust=1, vjust=1) + #geom_text
geom_text_repel(data=mnx[is.na(xm1),], aes(x=0.5*(lng+xd), y=0.5*(lat+yd), label=paste0(round(ddis,0), " km")), color = "black", size=3, hjust=1, vjust=1)
} else {
gf <- gf +
geom_segment(data=mnx, aes(x=lng, xend=xd, y=lat, yend=yd), size = 0.75, arrow = arrow(length = unit(0.08, "cm")), color="purple") +
geom_text_repel(data=mnx, aes(x=0.5*(lng+xd), y=0.5*(lat+yd), label=paste0(round(ddis,0), " km")), color = "black", size=3, hjust=1, vjust=1)
}
gf 